A Counterfactual Simulation Model of Causal Judgment
Tobias Gerstenberg
September 11, 2020
1 Setup
1.1 Load packages
library("knitr") # for knitting stuff
library("kableExtra") # for markdown tables
library("lme4") # for linear mixed effects models
library("broom.mixed") # for tidying mixed models results
library("brms") # Bayesian regression with Stan
library("corrr") # for tidy correlation matrix
library("xtable") # for latex tables
library("jsonlite") # for reading json files
library("janitor") # cleans stuff
library("Hmisc") # bootstrapped confidence intervals
library("tidybayes") # for Bayesian data analysis
library("png") # adding pngs to images
library("grid") # functions for dealing with images
library("plotly") # 3D scatter plot
library("egg") # for geom_custom
library("clValid") # for validating clustering solutions
# library("ggtern") # for ternary plots
library("patchwork") # for figure panels
library("tidyverse") # for wrangling, plotting, etc. 1.2 Helper functions
# setting some knitr options
opts_chunk$set(comment = "",
fig.show = "hold")
# set the ggplot theme
theme_set(theme_classic())
# suppress summarise() grouping warning
options(dplyr.summarise.inform = F)
# function for printing out html or latex tables
print_table = function(data, format = "html", digits = 2){
if(format == "html"){
data %>%
kable(digits = digits) %>%
kable_styling()
}else if(format == "latex"){
data %>%
xtable(digits = digits) %>%
print(include.rownames = F,
booktabs = T)
}
}
# root mean squared error
rmse = function(x, y){
return(sqrt(mean((x - y)^2)))
}
actual_counterfactual_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = clipindex,
y = value,
fill = colorindex,
group = model)) +
stat_summary(geom = "bar",
fun = "mean",
color = "black",
position = position_dodge(0.7),
width = 0.7) +
stat_summary(geom = "linerange",
fun.data = "mean_cl_boot",
position = position_dodge(0.7),
size = 1) +
geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
shape = 21,
color = "black",
size = 1,
alpha = 0.2) +
scale_fill_manual(values = c("red2", "green2", "black")) +
facet_grid(index_actual ~ index_counterfactual) +
geom_text(data = df.text %>%
drop_na(label),
mapping = aes(x = x, y = y, label = as.character(label)),
size = 6,
position = position_dodge(width = .7),
hjust = 0.5) +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 105),
expand = F,
clip = "off") +
labs(y = ylabel) +
theme(text = element_text(size = 20),
legend.position = "none",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
panel.spacing.y = unit(0.8, "cm"),
plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
axis.title.x = element_blank(),
panel.grid = element_blank(),
strip.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA))
# print(p)
}
actual_counterfactual_threeballs_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = ball,
y = value,
group = model,
fill = colorindex)) +
stat_summary(fun = "mean",
geom = "bar",
color = "black",
position = position_dodge(0.9),
width = 0.9) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
color = "black",
position = position_dodge(0.9)) +
geom_point(data = x %>%
mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(jitter.width = 0.3,
jitter.height = 0,
dodge.width = 0.9),
shape = 21,
color = "black",
size = 1,
alpha = 0.15) +
facet_wrap(~clip, ncol = 8) +
geom_text(data = df.text,
mapping = aes(x = x,
y = y,
label = as.character(label)),
size = 5,
position = position_dodge(width = .9),
hjust = 0.5,
fontface = "bold") +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 120),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25)) +
scale_fill_manual(values = c("red2", "green2", "black")) +
scale_color_manual(values = c("red2", "green2", "black")) +
labs(y = ylabel) +
theme(text = element_text(size = 16),
legend.position = "none",
axis.title.x = element_blank(),
panel.spacing.y = unit(0.3, "cm"),
panel.grid = element_blank(),
strip.background = element_blank(),
strip.text = element_blank(),
panel.border = element_rect(colour = "black", fill = NA))
}2 Experiment 1: 2 balls
2.1 Clips
2.2 Read in data
# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv",
header = T,
stringsAsFactors = F)
# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
header = T,
stringsAsFactors = F)
# Information about each clip
df.exp1.clipinfo = tibble(clip = 1:18,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 1),
outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1,
0, 1, 1, 0, 0, 1, 0, 1, 1),
index_actual = rep(c("actual miss",
"actual close",
"actual hit"), each = 6),
index_counterfactual = rep(rep(c("counterfactual miss",
"counterfactual close",
"counterfactual hit"),
each = 2),
3)) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp1.causal.long = df.exp1.causal$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:18, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
arrange(participant, clip, rating) %>%
left_join(df.exp1.clipinfo,
by = "clip")
df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:18, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
arrange(participant, clip, rating) %>%
left_join(df.exp1.clipinfo,
by = "clip")2.2.1 Demographics
# counterfactual condition
df.exp1.counterfactual %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 32 | 8.8 | 19 | 7.5 | 3.12 |
# causal condition
df.exp1.causal %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 33 | 10.3 | 21 | 9.65 | 6.58 |
2.3 Read in model predictions
2.3.1 Counterfactual condition
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]
df.exp1.model = files %>%
map_dfr(~ read.csv(str_c(path, .))) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp1.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup() %>%
left_join(df.exp1.clipinfo, by = "clip")
# find noisy simulation model that best predicts the mean counterfactual judgments by
# minimizing the sum of squared errors
df.exp1.counterfactual.model = df.exp1.model %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 -
df.exp1.counterfactual.means$rating_mean) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)2.3.2 Causal condition
df.exp1.causal.model = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(empirical = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.counterfactual.model %>%
select(-c(noise, sse)),
by = "clip") %>%
rename(simulation = prediction) %>%
mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))2.4 Counterfactual condition
2.4.1 Plots
set.seed(1)
df.plot = df.exp1.counterfactual.long %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.)/2)) %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss", "actual close", "actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")),
model = factor(model, levels = c("rating", "prediction"))) %>%
mutate(clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(x = df.plot, ylabel = "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
# width = 8,
# height = 6)
2.4.2 Stats
2.4.2.1 Model evaluation
df.exp1.counterfactual.means %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(noise = unique(noise),
simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
print_table()| noise | simulation_r | simulation_rmse |
|---|---|---|
| 1 | 0.97 | 10.95 |
2.5 Causal condition
2.5.1 Plots
set.seed(1)
model.name = c("empirical")
# model.name = c("simulation")
df.plot = df.exp1.causal.long %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, simulation, empirical),
names_to = "model",
values_to = "value") %>%
filter(model %in% c(model.name, "rating")) %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(
colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index.actual = factor(index_actual,
levels = c("actual miss", "actual close", "actual hit")),
index.counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit"))) %>%
mutate(clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)) # swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(df.plot, "causal judgment")
print(p)
# ggsave("../../figures/plots/exp1_causal_bars.pdf",
# width = 8,
# height = 6)
2.5.2 Stats
2.5.2.1 Model evaluation
df.exp1.causal.long %>%
mutate(rating = abs(rating)) %>%
group_by(clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(empirical_r = cor(rating, empirical),
empirical_rmse = rmse(rating, empirical),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
print_table()| empirical_r | empirical_rmse | simulation_r | simulation_rmse |
|---|---|---|---|
| 0.96 | 8.57 | 0.93 | 15.15 |
2.5.2.2 Bayesian regression
2.5.2.3 Specific hypotheses
df.exp1.causal.posterior = df.exp1.causal.long %>%
group_by(clip, index_actual, index_counterfactual, outcome_actual) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp1_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, .value, .draw) %>%
pivot_wider(names_from = clip,
values_from = .value)
func_posterior_difference = function(data, clip1, clip2){
data %>%
mutate(difference = .data[[clip1]] - .data[[clip2]]) %>%
pull(difference) %>%
mean_hdci() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>%
print()
}
func_posterior_difference(data = df.exp1.causal.posterior,
clip1 = "8",
clip2 = "13") difference
1 (0.26 [-6.51, 7.31])2.5.3 Tables
fit.brm_exp1_causal %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
# print_table(format = "latex")
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | -18.33 | 3.36 | -24.86 | -11.66 |
| fixed | cond | index_actualactualclose | 4.31 | 3.40 | -2.54 | 10.94 |
| fixed | cond | index_actualactualhit | 4.06 | 4.84 | -5.51 | 13.27 |
| fixed | cond | index_counterfactualcounterfactualclose | -17.52 | 3.37 | -24.36 | -10.84 |
| fixed | cond | index_counterfactualcounterfactualhit | -61.37 | 4.38 | -69.95 | -52.89 |
| fixed | cond | outcome_actual | 92.36 | 5.11 | 82.21 | 102.38 |
3 Experiment 2: Brick & Teleport
3.1 Clips
3.2 Read in data
df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
header = T,
stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
header = T,
stringsAsFactors = F)
# Information about each clip
df.exp2.clipinfo = tibble(clip = 1:20,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0,
1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
index_actual = c(rep(c("actual miss",
"actual close",
"actual hit"),
each = 6),
rep("practice", 2)),
index_counterfactual = c(rep(rep(c("counterfactual miss",
"counterfactual close",
"counterfactual hit"),
each = 2),
3),
rep("practice", 2))) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp2.causal_first.long = df.exp2.causal_first$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:40, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(question = rep(rep(c("causal", "counterfactual"), each = 20),
max(participant)),
condition = "causal_first") %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:40, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(question = rep(rep(c("counterfactual", "causal"), each = 20),
max(participant)),
condition = "counterfactual_first") %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
bind_rows(df.exp2.counterfactual_first.long %>%
mutate(participant = participant +
max(df.exp2.causal_first.long$participant)))3.2.1 Demographics
# counterfactual condition
df.exp2.counterfactual_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time) %>%
bind_rows(df.exp2.causal_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time)) %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 82 | 35 | 12.2 | 45 | 21.25 | 5.11 |
3.3 Read in model predictions
3.3.1 Counterfactual condition
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]
df.exp2.model = files %>%
map_dfr(~ read.csv(str_c(path, .))) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp2.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp2.counterfactual.means = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
# mutate(rating = abs(rating)) %>%
group_by(clip) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.clipinfo, by = "clip")
# find noisy simulation model that best predicts the mean counterfactual judgments by
# minimizing the sum of squared errors
df.exp2.counterfactual.model = df.exp2.model %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 - df.exp2.counterfactual.means$rating) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)3.3.2 Causal condition
# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
group_by(clip, condition) %>%
summarize(rating = mean(rating)) %>%
pivot_wider(names_from = condition,
values_from = rating) %>%
left_join(df.exp2.combined.long %>%
filter(question == "counterfactual", clip <= 18) %>%
group_by(clip) %>%
summarize(combined = mean(rating)),
by = "clip") %>%
left_join(df.exp2.counterfactual.model,
by = "clip") %>%
select(-c(sse, noise)) %>%
rename(simulation = prediction) %>%
mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined),
.funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))3.4 Counterfactual condition
3.4.1 Plots
set.seed(1)
df.plot = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")),
model = factor(model,
levels = c("rating", "prediction")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"),
each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
# width = 8,
# height = 6)
3.4.2 Stats
3.4.2.1 Model evaluation
df.exp2.counterfactual.means %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(noise = unique(noise),
simulation_r = cor(rating, prediction),
simulation_rmse = rmse(rating, prediction)) %>%
print_table()| noise | simulation_r | simulation_rmse |
|---|---|---|
| 0.9 | 0.96 | 13.46 |
3.4.2.2 Bayesian regressions
3.4.2.2.1 Check for potential order effects
# random intercepts & slopes model
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition *
(index_actual + index_counterfactual) +
(1 + index_actual +
index_counterfactual | participant),
data = df.exp2.combined.long %>%
na.omit() %>%
filter(question == "counterfactual"),
cores = 2,
seed = 1,
file = "cache/brm_exp2_counterfactual")3.4.3 Tables
fit.brm_exp2_counterfactual %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
# print_table(format = "latex")
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | 13.38 | 2.58 | 8.30 | 18.38 |
| fixed | cond | conditioncounterfactual_first | -1.11 | 3.65 | -7.99 | 6.04 |
| fixed | cond | index_actualactualclose | -1.48 | 2.78 | -6.79 | 3.95 |
| fixed | cond | index_actualactualhit | -4.42 | 2.40 | -9.15 | 0.20 |
| fixed | cond | index_counterfactualcounterfactualclose | 37.94 | 3.33 | 31.43 | 44.60 |
| fixed | cond | index_counterfactualcounterfactualhit | 73.80 | 4.17 | 65.54 | 81.99 |
| fixed | cond | conditioncounterfactual_first:index_actualactualclose | -0.55 | 3.84 | -8.29 | 6.79 |
| fixed | cond | conditioncounterfactual_first:index_actualactualhit | -0.85 | 3.41 | -7.60 | 5.87 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualclose | 1.23 | 4.79 | -8.37 | 10.75 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualhit | -0.02 | 5.92 | -11.96 | 11.59 |
3.5 Causal condition
3.5.1 Plots
set.seed(1)
func_exp2_causal_plot = function(condition.name, model.name){
df.plot = df.exp2.combined.long %>%
filter(question == "causal", clip <= 18) %>%
filter(condition %in% condition.name) %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c("rating", all_of(model.name)),
names_to = "model",
values_to = "value") %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"),
each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
actual_counterfactual_plot(df.plot, "causal judgment")
}
plot.list = list(func_exp2_causal_plot(condition.name = "counterfactual_first",
model.name = "combined"),
func_exp2_causal_plot(condition.name = "causal_first",
model.name = "combined"))
wrap_plots(plot.list, ncol = 2)
# condition.name = "counterfactual_first"
# condition.name = "causal_first"
# model.name = "combined"
#
# func_exp2_causal_plot(condition.name = condition.name,
# model.name = model.name)
# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
# width = 8,
# height = 6)
Figure 3.1: Left panel: Counterfactual judgments first; right panel: Causal judgments first
3.5.2 Stats
3.5.2.1 Model evaluation
df.data = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
mutate(rating = abs(rating)) %>%
group_by(condition,
clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual"))
df.data %>%
group_by(condition) %>%
summarize(empirical_r = cor(rating, combined),
empirical_rmse = rmse(rating, combined)) %>%
print_table()| condition | empirical_r | empirical_rmse |
|---|---|---|
| causal_first | 0.92 | 13.78 |
| counterfactual_first | 0.99 | 5.27 |
df.data %>%
summarize(empirical_r = cor(rating, combined),
empirical_rmse = rmse(rating, combined),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
print_table()| empirical_r | empirical_rmse | simulation_r | simulation_rmse |
|---|---|---|---|
| 0.95 | 10.43 | 0.92 | 19.09 |
3.5.2.2 Bayesian regression
# random intercepts + slopes model
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition *
(index_actual + index_counterfactual + outcome_actual) +
(1 + index_actual + index_counterfactual +
outcome_actual | participant),
data = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18),
cores = 2,
seed = 1,
file = "cache/brm_exp2_causal")3.5.2.3 Specific hypotheses
df.exp2.causal.posterior = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
group_by(clip, index_actual, index_counterfactual,
outcome_actual, condition) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp2_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, condition, .value, .draw) %>%
pivot_wider(names_from = clip,
values_from = .value)
# difference between clips 13 and 14 versus 8
df.exp2.causal.posterior %>%
mutate(difference = (`13` + `14`) / 2 - `8`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) difference
1 (-0.55 [-12.09, 10.97])# difference between clips 5 and 6 versus 11
df.exp2.causal.posterior %>%
mutate(difference = (`5` + `6`) / 2 - `11`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) difference
1 (-5.48 [-13.8, 2.73])3.5.3 Tables
3.5.3.1 Bayesian regression
fit.brm_exp2_causal %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
# print_table(format = "latex")
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | -24.07 | 3.50 | -30.87 | -17.24 |
| fixed | cond | conditioncounterfactual_first | 11.47 | 4.91 | 1.62 | 21.22 |
| fixed | cond | index_actualactualclose | 8.17 | 3.44 | 1.47 | 14.79 |
| fixed | cond | index_actualactualhit | 13.42 | 4.90 | 3.74 | 22.69 |
| fixed | cond | index_counterfactualcounterfactualclose | -30.20 | 3.72 | -37.61 | -23.19 |
| fixed | cond | index_counterfactualcounterfactualhit | -50.29 | 4.60 | -59.57 | -41.47 |
| fixed | cond | outcome_actual | 98.53 | 5.40 | 87.81 | 108.98 |
| fixed | cond | conditioncounterfactual_first:index_actualactualclose | -5.39 | 4.94 | -15.00 | 4.35 |
| fixed | cond | conditioncounterfactual_first:index_actualactualhit | -16.97 | 7.02 | -30.40 | -3.33 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualclose | -8.13 | 5.17 | -18.18 | 2.17 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualhit | -25.51 | 6.44 | -38.16 | -12.98 |
| fixed | cond | conditioncounterfactual_first:outcome_actual | 10.72 | 7.60 | -4.16 | 26.02 |
4 Experiment 3: 3 balls
4.1 Clips
4.2 Read in data
# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1,
0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
clipindex = rep(1:2, 16))
# COUNTERFACTUAL JUDGMENTS
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
header = T,
stringsAsFactors = F)
# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("gender", "age", "difficulty",
"smoothness", "time", "condition"),
n()/6),
participant = rep(1:(n()/6), each = 6)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
select(participant, condition, gender, age, time, smoothness, difficulty)
# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("clip", "rating"), n()/2),
participant = rep(1:(n()/64), each = 64),
order = rep(rep(1:32, each = 2), n()/64)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
left_join(df.exp3.counterfactual.demographics %>%
select(participant, ball = condition),
by = "participant") %>%
select(participant, ball, clip, everything()) %>%
arrange(participant, clip)
# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
header = T,
stringsAsFactors = F)
# demographics
df.exp3.causal.demographics = df.exp3.causal$extra %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("gender", "age", "difficulty", "smoothness", "time",
"counterbalance", "replayType", "experiment"),
n()/8),
participant = rep(1:(n()/8), each = 8)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
select(participant, gender, age, counterbalance, time, smoothness, difficulty)
# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
participant = rep(1:(n()/(32*6)), each = (32*6)),
order = rep(rep(1:32, each = 6), n()/(32*6))) %>%
filter(!name %in% c("x", "y", "z")) %>%
pivot_wider(names_from = name,
values_from = value) %>%
left_join(df.exp3.causal.demographics %>%
select(participant, counterbalance),
by = "participant") %>%
mutate(A = ifelse(counterbalance == 1, rating1, rating2),
B = ifelse(counterbalance == 1, rating2, rating1)) %>%
select(-rating1, -rating2) %>%
pivot_longer(cols = c(A, B),
names_to = "ball",
values_to = "rating") %>%
select(participant, clip, ball, rating, order) %>%
arrange(participant, clip, ball)4.2.1 Demographics
# counterfactual condition
df.exp3.counterfactual.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 80 | 33 | 10.1 | 34 | 18.08 | 4.63 |
# causal condition
df.exp3.causal.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 34 | 10.5 | 21 | 21.19 | 4.96 |
4.3 Read in model predictions
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]
df.exp3.model = files %>%
map_dfr(~ read_csv(str_c(path, .))) %>%
rename(clip = trial)4.3.1 Counterfactual condition
# calculate mean counterfactual judgments
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()
# find noisy simulation model that best predicts the mean counterfactual
# judgments by minimizing the sum of squared errors
df.exp3.model.counterfactual = df.exp3.model %>%
group_by(noise) %>%
select(clip, contains("whether"), noise) %>%
pivot_longer(cols = c(A_whether, B_whether),
names_to = "ball",
values_to = "prediction") %>%
mutate(ball = str_remove(ball, "_whether")) %>%
arrange(noise, clip, ball) %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 -
df.exp3.counterfactual.means$rating_mean) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)4.3.2 Causal condition
# calculate mean causal judgments
df.exp3.causal.means = df.exp3.causal.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()
df.exp3.model.causal = df.exp3.model %>%
# take into account difference-making
mutate_at(.vars = vars(A_how:A_robust),
.funs = list(~ . * A_difference)) %>%
mutate_at(.vars = vars(B_how:B_robust),
.funs = list(~ . * B_difference)) %>%
# choose model based on best fit with counterfactual data
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
pivot_longer(cols = A_difference:B_robust,
names_to = c("ball", "aspect"),
names_sep = "_",
values_to = "value") %>%
pivot_wider(names_from = aspect,
values_from = value) %>%
select(clip, ball:robust)4.3.3 Heuristic model
l.features = fromJSON("data/features.json")
df.features.initial = l.features[["appearance"]] %>%
cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
"ballA_velx",
"ballA_vely",
"ballB_velx",
"ballB_vely",
"ballE_velx",
"ballE_vely")) %>%
mutate(clip = 1:nrow(.)) %>%
pivot_longer(cols = starts_with("ball")) %>%
separate(name, into = c("ball", "property")) %>%
pivot_wider(names_from = property, values_from = value) %>%
mutate(ball = str_remove_all(ball, pattern = "ball"))
# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")
for (i in 1:32) {
df.tmp = data.frame()
for (j in 1:length(ballnames)) {
ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
if (ncollisions > 0) {
tmp = data.frame(
ball = ballnames[j],
object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
ncollision = 1:ncollisions
)
df.tmp = df.tmp %>%
rbind(tmp)
}
}
df.features.collisions = df.features.collisions %>%
rbind(df.tmp %>%
mutate(clip = i) %>%
select(clip, ball, everything()))
}
# find end of clip
tmp = df.features.collisions %>%
group_by(clip) %>%
filter(ball == "ballE",
str_detect(object, "goal|Left|Right")) %>%
group_by(clip) %>%
mutate(endclip = max(time)) %>%
select(clip, endclip) %>%
ungroup() %>%
distinct()
# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
left_join(tmp, by = "clip") %>%
mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
group_by(clip) %>%
filter(time <= endclip)
# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
filter(ball == "E") %>%
mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
select(clip, E.moving)
# contact with ball E
tmp.contactE = df.features.collisions %>%
filter(object == "ballE") %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
mutate(contactE = 1) %>%
select(clip, ball, contactE) %>%
group_by(clip, ball) %>%
summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate(contactE = ifelse(is.na(contactE), 0, contactE))
# number of collisions
tmp.ncollisions = df.features.collisions %>%
filter(str_detect(object, "ball"),
ball != "ballE") %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
group_by(clip, ball) %>%
count() %>%
rename(ncollision = n)
# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
filter(str_detect(object, "ball"),
ball != "ballE") %>%
count(clip, ball, object) %>%
mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
group_by(clip) %>%
summarize(AE = any(AE %>% as.logical()) * 1,
BE = any(BE %>% as.logical()) * 1) %>%
mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
select(clip, A, B) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "exclusive")
# impact
func_angle = function(x, y) {
dot.prod = x %*% y
norm.x = norm(x, type = "2")
norm.y = norm(y, type = "2")
theta = acos(dot.prod / (norm.x * norm.y))
as.numeric(theta)
}
# impact on ball E
tmp.impactE = df.features.collisions %>%
rowwise() %>%
filter(str_detect(object, "ball"),
ball == "ballE") %>%
mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely),
c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff),
abs(atan2(post.vely - pre.vely,
post.velx - pre.velx)),
direction.diff)) %>%
ungroup() %>%
group_by(clip, object) %>%
summarize(speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)) %>%
rename(ball = object) %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate_at(.vars = vars(speed.diff, direction.diff),
.funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(E.speed.diff = speed.diff,
E.direction.diff = direction.diff)
# total impact
tmp.impactTotal = df.features.collisions %>%
rowwise() %>%
filter(str_detect(object, "ballA|ballB")) %>%
mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely),
c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff),
abs(atan2(post.vely - pre.vely,
post.velx - pre.velx)),
direction.diff)) %>%
ungroup() %>%
group_by(clip, object) %>%
summarize(speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)) %>%
rename(ball = object) %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate_at(.vars = vars(speed.diff, direction.diff),
.funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(total.speed.diff = speed.diff,
total.direction.diff = direction.diff)
# transfer of force
tmp.transfer = df.features.collisions %>%
filter(str_detect(ball, "ballA|ballB"),
str_detect(object, "ball")) %>%
mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
arrange(clip, time) %>%
group_by(clip) %>%
summarize(A = any(!is.na(AE)),
B = any(!is.na(BE)),
A = ifelse(any(!is.na(AB)) &
any(!is.na(BE)) &
max(AB, na.rm = T) < min(BE, na.rm = T),
T,
A),
B = ifelse(any(!is.na(AB)) &
any(!is.na(AE)) &
max(AB, na.rm = T) < min(AE, na.rm = T),
T,
B)) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "transfer") %>%
arrange(clip, ball)
# collect features
df.features = df.features.initial %>%
filter(ball != "E") %>%
mutate(ball = factor(ball)) %>%
mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
speed = abs(velx) + abs(vely)) %>%
left_join(tmp.contactE) %>%
left_join(tmp.impactE) %>%
left_join(tmp.impactTotal) %>%
left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
left_join(tmp.movingE) %>%
left_join(tmp.exclusive) %>%
select(-c(appearance, velx, vely))
# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])4.4 Counterfactual condition
4.4.1 Plots
4.4.1.1 Bar graph
set.seed(1)
df.plot = df.exp3.counterfactual.long %>%
left_join(df.exp3.model.counterfactual %>%
select(clip, ball, prediction) %>%
mutate(ball = as.factor(ball)),
by = c("clip", "ball")) %>%
left_join(df.exp3.clipinfo,
by = "clip") %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = 1,
colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1,
2,
colorindex),
colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1,
2,
colorindex),
colorindex = ifelse(model == "prediction", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "prediction"))) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.5,
colorindex = NA)
p = actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
# width = 12,
# height = 6)
4.4.2 Stats
4.4.2.1 Model evaluation
# best fitting model
df.exp3.counterfactual.means %>%
left_join(df.exp3.model.counterfactual,
by = c("clip", "ball")) %>%
summarize(simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.86 | 19.84 |
# deterministic model
df.exp3.counterfactual.means %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_a, outcome_b) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "outcome") %>%
mutate(ball = factor(ball,
levels = c("outcome_a", "outcome_b"),
labels = c("B", "A"))),
by = c("clip", "ball")) %>%
mutate(outcome = outcome * 100) %>%
summarize(simulation_r = cor(rating_mean, outcome),
simulation_rmse = rmse(rating_mean, outcome)) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.83 | 30.28 |
4.5 Causal condition
4.5.1 Stats
4.5.1.1 Bayesian mixed effects model
df.data = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
# using whether-causation from the model
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether +
(1 + whether | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_w")
fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how +
(1 + whether + how | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_wh")
fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient +
(1 + whether + how + sufficient | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_whs")4.5.1.2 Model comparison
fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object. elpd_diff se_diff
fit.brm_exp3_causal_whs 0.0 0.0
fit.brm_exp3_causal_wh -125.2 15.6
fit.brm_exp3_causal_w -426.2 31.1 4.5.1.3 Heuristic model
- z-scored predictors and invidiual responses
# Regression based on features
df.regression.features = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
left_join(df.features %>%
mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
by = c("clip", "ball")) %>%
clean_names() %>%
mutate(e_moving = 1 - e_moving)
# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving +
speed +
contact_e +
e_speed_diff +
e_direction_diff +
total_speed_diff +
total_direction_diff +
transfer +
e_moving +
exclusive + (1 | participant),
prior = priors,
data = df.regression.features %>%
select(-c(clip, ball, order, clipindex,
contains("outcome"))),
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant)
Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 8.53 1.22 6.44 11.18 1.00 1163 1936
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 49.79 1.46 46.91 52.74 1.00 1129
moving 0.22 0.22 0.01 0.81 1.00 3399
speed 2.06 0.85 0.44 3.67 1.00 2319
contact_e 0.38 0.35 0.01 1.29 1.00 3201
e_speed_diff 0.12 0.12 0.00 0.46 1.00 4435
e_direction_diff 1.04 0.74 0.04 2.71 1.00 2202
total_speed_diff 2.21 0.95 0.41 4.17 1.00 2208
total_direction_diff 3.99 0.91 2.14 5.68 1.00 2718
transfer 15.59 0.79 14.06 17.14 1.00 3555
e_moving 0.18 0.17 0.01 0.62 1.00 3296
exclusive 4.39 0.73 2.97 5.81 1.00 3526
Tail_ESS
Intercept 1854
moving 1547
speed 1266
contact_e 2020
e_speed_diff 2268
e_direction_diff 1685
total_speed_diff 1393
total_direction_diff 1742
transfer 2992
e_moving 1738
exclusive 2538
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 33.21 0.46 32.33 34.10 1.00 5313 2845
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).4.5.1.4 Predictions for mean causal judgments
func_fit_data = function(fit){
fit %>%
fitted(newdata = df.exp3.model.causal %>%
left_join(df.features %>%
mutate_at(.vars = vars(moving:exclusive),
.funs = ~ scale(.)[,]),
by = c("clip", "ball")) %>%
clean_names() %>%
mutate(e_moving = 1 - e_moving),
re_formula = NA) %>%
as_tibble() %>%
pull(Estimate)
}
df.exp3.causal.regression = df.exp3.causal.means %>%
mutate(w = func_fit_data(fit.brm_exp3_causal_w),
wh = func_fit_data(fit.brm_exp3_causal_wh),
whs = func_fit_data(fit.brm_exp3_causal_whs),
heuristic = func_fit_data(fit.brm_exp3_causal_heuristic))4.5.1.5 Model evaluation
prediction = "whs"
df.exp3.causal.regression %>%
summarize(simulation_r = cor(rating_mean, .[[prediction]]),
simulation_rmse = rmse(rating_mean, .[[prediction]])) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.87 | 13.07 |
4.5.1.6 Individual participant regressions
df.fit = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
# priors
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline =
brm(formula = rating ~ 1,
data = df.fit %>%
filter(participant == 1),
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_baseline"))
fit.brm_exp3_causal_individual_w =
brm(formula = rating ~ 1 + whether,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_w"))
fit.brm_exp3_causal_individual_wh =
brm(formula = rating ~ 1 + whether + how,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_wh"))
fit.brm_exp3_causal_individual_whs =
brm(formula = rating ~ 1 + whether + how + sufficient,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_whs"))
# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
group_by(participant) %>%
nest() %>%
ungroup() %>%
# fit model for each participant
mutate(fit_baseline = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_baseline,
newdata = .x)),
fit_w = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_w,
newdata = .x)),
fit_wh = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_wh,
newdata = .x)),
fit_whs = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_whs,
newdata = .x))) %>%
mutate(fit_baseline = map(.x = fit_baseline,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_w = map(.x = fit_w,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_wh = map(.x = fit_wh,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_whs = map(.x = fit_whs,
~ add_criterion(.x, criterion = c("loo", "waic"))),
r_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
model_comparison = pmap(.l = list(baseline = fit_baseline,
w = fit_w,
wh = fit_wh,
whs = fit_whs),
.f = ~ loo_compare(..1, ..2, ..3, ..4)),
best_model = map_chr(.x = model_comparison,
.f = ~ rownames(.) %>% .[1]),
best_model = factor(best_model,
levels = c("..1", "..2", "..3", "..4"),
labels = c("baseline", "w", "wh", "whs")))
save(list = c("df.exp3.causal.individual_fit"),
file = "data/exp3_causal_individual_fit.RData")4.5.1.6.1 Load individual participant regressions
load("data/exp3_causal_individual_fit.RData")
# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>%
count(best_model) %>%
print_table()| best_model | n |
|---|---|
| wh | 2 |
| whs | 39 |
4.5.1.6.2 Model performance on individual participant level
df.exp3.causal.individual_fit %>%
select(participant, contains("r_"), contains("rmse_")) %>%
pivot_longer(cols = -participant,
names_to = c("measure", "model"),
names_sep = "_",
values_to = "fit") %>%
group_by(model, measure) %>%
summarize(quantiles = quantile(fit, probs = c(0.05, 0.5, 0.95)),
prob = c(0.05, 0.5, 0.95)) %>%
pivot_wider(names_from = prob,
values_from = quantiles) %>%
mutate(across(.fns = ~ round(., 2))) %>%
print_table()| model | measure | 0.05 | 0.5 | 0.95 |
|---|---|---|---|---|
| w | r | 0.18 | 0.43 | 0.57 |
| w | rmse | 29.33 | 36.03 | 43.62 |
| wh | r | 0.35 | 0.60 | 0.78 |
| wh | rmse | 23.82 | 32.19 | 40.60 |
| whs | r | 0.39 | 0.67 | 0.79 |
| whs | rmse | 22.69 | 31.17 | 38.12 |
4.5.1.6.3 Clusters of individual participants’ responses
set.seed(1)
df.cluster_whs = fit.brm_exp3_causal_whs %>%
ranef() %>%
.$participant %>%
as_tibble() %>%
select(contains("Estimate")) %>%
set_names(tolower(str_replace_all(names(.), "Estimate.", ""))) %>%
mutate(participant = 1:n()) %>%
relocate(participant)
df.cluster_whs = df.cluster_whs %>%
mutate(cluster = df.cluster_whs %>%
select(-participant) %>%
kmeans(centers = 2) %>%
.$cluster)
df.cluster_whs %>%
print_table()| participant | intercept | whether | how | sufficient | cluster |
|---|---|---|---|---|---|
| 1 | 8.30 | 0.44 | -5.36 | -0.84 | 1 |
| 2 | 0.03 | -20.89 | 32.67 | -20.83 | 2 |
| 3 | 2.27 | -8.27 | 26.27 | -9.42 | 2 |
| 4 | -0.66 | 7.98 | -2.96 | 10.74 | 1 |
| 5 | -0.75 | 2.47 | -5.92 | 0.20 | 1 |
| 6 | 5.09 | -15.71 | 22.42 | -18.09 | 2 |
| 7 | -0.17 | 1.71 | -3.13 | 0.22 | 1 |
| 8 | -3.55 | -2.49 | 4.88 | -0.21 | 2 |
| 9 | -4.15 | -5.31 | 9.38 | -2.57 | 2 |
| 10 | 2.10 | -9.99 | 3.25 | -10.96 | 2 |
| 11 | 2.10 | 7.99 | -10.81 | 5.96 | 1 |
| 12 | -2.00 | 15.83 | -25.89 | 18.27 | 1 |
| 13 | 3.03 | 4.98 | 0.72 | 2.29 | 1 |
| 14 | 3.17 | 4.02 | -11.54 | 1.56 | 1 |
| 15 | -5.07 | -13.56 | -3.49 | -11.90 | 2 |
| 16 | -0.89 | 23.98 | -17.31 | 25.30 | 1 |
| 17 | 0.35 | -4.04 | 2.84 | -3.79 | 2 |
| 18 | 3.05 | -0.15 | -10.15 | -2.03 | 1 |
| 19 | -2.30 | -7.53 | 4.17 | -7.76 | 2 |
| 20 | 8.44 | 20.15 | -11.67 | 17.52 | 1 |
| 21 | -1.69 | 1.54 | -12.89 | 2.37 | 1 |
| 22 | -5.38 | 1.33 | -3.82 | 2.56 | 1 |
| 23 | -0.53 | 1.24 | -7.19 | 3.29 | 1 |
| 24 | -0.43 | 1.91 | -0.12 | 2.55 | 1 |
| 25 | 2.88 | 12.33 | -12.26 | 9.65 | 1 |
| 26 | 0.46 | 12.55 | 1.40 | 11.74 | 1 |
| 27 | -4.46 | 16.00 | -20.07 | 21.43 | 1 |
| 28 | 2.21 | -8.24 | -1.57 | -10.35 | 2 |
| 29 | -0.78 | 7.61 | -4.96 | 7.93 | 1 |
| 30 | -3.78 | -21.34 | 18.53 | -19.98 | 2 |
| 31 | -6.59 | 7.50 | -12.71 | 12.17 | 1 |
| 32 | -1.70 | 2.50 | 3.37 | 6.52 | 1 |
| 33 | 3.29 | 13.81 | -10.56 | 13.26 | 1 |
| 34 | -4.20 | -12.18 | 9.83 | -10.61 | 2 |
| 35 | 0.41 | -16.26 | 0.37 | -18.28 | 2 |
| 36 | -1.55 | -27.69 | 47.62 | -28.75 | 2 |
| 37 | 10.64 | 13.67 | 0.43 | 8.39 | 1 |
| 38 | 1.35 | 9.88 | -7.13 | 9.38 | 1 |
| 39 | -5.71 | -10.20 | 7.58 | -10.99 | 2 |
| 40 | -6.18 | -7.24 | 8.60 | -6.13 | 2 |
| 41 | 1.97 | -2.12 | -1.61 | -0.67 | 1 |
# A tibble: 2 x 2
cluster n
<int> <int>
1 1 25
2 2 164.5.1.6.3.1 Check that clustering is stable
A 2-cluster solution yields a good result according to a number of different validation measures (see here for more details on the different measures).
tmp = df.cluster_whs %>%
select(-participant, -cluster)
fit = clValid(obj = as.matrix(tmp),
nClust = 2:10,
clMethods = c("kmeans"),
validation = c("internal", "stability"))Warning in clValid(obj = as.matrix(tmp), nClust = 2:10, clMethods =
c("kmeans"), : rownames for data not specified, using 1:nrow(data)
Clustering Methods:
kmeans
Cluster sizes:
2 3 4 5 6 7 8 9 10
Validation Measures:
2 3 4 5 6 7 8 9 10
kmeans APN 0.0991 0.1419 0.1438 0.1181 0.2214 0.1734 0.2094 0.1964 0.2269
AD 20.3268 15.8118 13.3966 11.5686 11.4495 10.0034 9.6513 9.2184 8.8807
ADM 7.3494 4.1130 4.3773 2.2086 3.7885 3.1164 3.5423 3.9327 3.9086
FOM 8.3657 7.0101 6.0226 5.4730 5.1633 4.9791 4.8082 4.8092 4.9224
Connectivity 4.4262 15.5893 20.1163 22.5329 25.0329 28.1413 34.4710 42.6353 44.9687
Dunn 0.1046 0.1626 0.1687 0.2509 0.2509 0.2550 0.1815 0.2524 0.2524
Silhouette 0.4691 0.3871 0.3868 0.3738 0.3522 0.3255 0.2988 0.2846 0.2774
Optimal Scores:
Score Method Clusters
APN 0.0991 kmeans 2
AD 8.8807 kmeans 10
ADM 2.2086 kmeans 5
FOM 4.8082 kmeans 8
Connectivity 4.4262 kmeans 2
Dunn 0.2550 kmeans 7
Silhouette 0.4691 kmeans 2 4.5.2 Plots
4.5.2.1 Bar graph
set.seed(1)
model_index = "whs"
# model predictions
model_prediction = list(fit.brm_exp3_causal_w,
fit.brm_exp3_causal_wh,
fit.brm_exp3_causal_whs) %>%
map_dfr(~ fitted(., df.exp3.causal.means %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")),
re_formula = NA) %>%
as_tibble()) %>%
mutate(ball = rep(c("A", "B"), n()/2),
clip = rep(rep(1:32, each = 2), 3),
model = rep(c("w", "wh", "whs"),
each = 64))
df.plot = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_both),
by = c("clip")) %>%
left_join(model_prediction %>%
filter(model == model_index) %>%
select(model = Estimate, clip, ball),
by = c("clip", "ball")) %>%
pivot_longer(cols = c(rating, model),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = 1,
colorindex = ifelse(model == "rating" &
outcome_both == 0, 1, colorindex),
colorindex = ifelse(model == "rating" &
outcome_both == 1, 2, colorindex),
colorindex = ifelse(model == "model", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "model"))) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.5,
colorindex = NA)
p = actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
print(p)
# ggsave("../../figures/plots/exp3_causal_bars.pdf",
# width = 12,
# height = 6)
4.5.2.2 Selection
check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(3, 7, 15, 16, 23)) %>%
group_by(clip, ball) %>%
summarize(mean = mean(rating),
low = smean.cl.boot(rating)[2],
high = smean.cl.boot(rating)[3]) %>%
left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
by = c("clip", "ball")) %>%
ungroup() %>%
pivot_longer(cols = c(mean, w, wh, whs),
names_to = "index",
values_to = "value") %>%
mutate_at(.vars = vars(low, high),
.funs = ~ ifelse(index == "mean", ., NA)) %>%
mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain",
"double prevention",
"joint causation",
"overdetermination",
"preemption")))
df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
value = -10,
index = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
index = NA,
value = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = value,
group = index,
fill = index)) +
geom_bar(stat = "identity", color = "black",
position = position_dodge(0.9),
width = 0.9) +
geom_errorbar(mapping = aes(ymin = low, ymax = high),
width = 0,
size = 1,
position = position_dodge(0.9)) +
annotation_custom(grob = ball_a,
xmin = 0.5,
xmax = 1.5,
ymin = -25,
ymax = -7) +
annotation_custom(grob = ball_b,
xmin = 1.5,
xmax = 2.5,
ymin = -25,
ymax = -7) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_grid(~clip) +
scale_fill_grey(start = 1,
end = 0,
labels = c("mean rating ", expression(CSM[W] ~ " "),
expression(CSM[WH] ~ " "),
expression(CSM[WHS] ~ " "))) +
labs(y = "causal responsibility", fill = "") +
coord_cartesian(clip = "off") +
theme_bw() +
theme(legend.text.align = 0,
text = element_text(size = 20),
panel.grid = element_blank(),
legend.position = "bottom",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
legend.text = element_text(size = 30),
strip.text = element_text(size = 30),
legend.spacing = unit(0.5, "cm"),
legend.background = element_rect(fill = "transparent"),
legend.margin = margin(t = 5, unit = "cm"),
plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
print(p)
# ggsave("../../figures/plots/exp3_selection_bars.pdf",
# width = 20,
# height = 10,
# device = cairo_pdf)
4.5.2.3 Scatter plot
func_scatterplot = function(model){
if(model == "heuristic"){
xlabel = "Heuristic"
}else{
xlabel = bquote(CSM[.(toupper(model))])
}
l.models = list(w = fit.brm_exp3_causal_w,
wh = fit.brm_exp3_causal_wh,
whs = fit.brm_exp3_causal_whs,
heuristic = fit.brm_exp3_causal_heuristic)
df.data = df.exp3.causal.means %>%
left_join(df.exp3.model.causal, by = c("clip", "ball")) %>%
left_join(df.regression.features, by = c("clip", "ball"))
df.plot = l.models[[model]] %>%
fitted(newdata = df.data,
re_formula = NA) %>%
as_tibble() %>%
clean_names() %>%
bind_cols(df.data)
p = ggplot(data = df.plot,
mapping = aes(x = estimate,
y = rating_mean)) +
geom_abline(intercept = 0,
slope = 1,
linetype = 2) +
geom_smooth(aes(y = estimate,
ymin = q2_5,
ymax = q97_5),
stat = "identity",
color = "black") +
geom_linerange(size = 0.75,
mapping = aes(ymin = rating_low,
ymax = rating_high),
color = "gray50") +
geom_point(size = 2) +
scale_color_manual(values = c("black", "#e41a1c"),
guide = F) +
coord_fixed(xlim = c(0, 100),
ylim = c(0, 100)) +
labs(x = xlabel,
y = "responsibility judgment") +
annotate(geom = "text",
label = str_c(
"r = ", cor(df.plot$estimate, df.plot$rating_mean) %>%
round(2) %>%
as.character() %>%
str_sub(start = 2),
"\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>%
round(2)),
x = 5,
y = 95,
size = 6,
hjust = 0) +
scale_x_continuous(breaks = seq(0, 100, 25)) +
scale_y_continuous(breaks = seq(0, 100, 25)) +
theme(text = element_text(size = 20))
return(p)
}
# for creating and saving an individual scatter plot
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"
# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
# width = 5,
# height = 5)
list(func_scatterplot(model = "w"),
func_scatterplot(model = "wh"),
func_scatterplot(model = "whs"),
func_scatterplot(model = "heuristic")) %>%
wrap_plots(ncol = 2)
4.5.2.4 Individual participant variance for a selection of clips (clustered)
check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(7, 23, 3, 15, 16)) %>%
mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain", "double prevention", "joint causation",
"overdetermination", "preemption")))
df.plot = df.plot %>%
left_join(df.cluster_whs %>%
group_by(cluster) %>%
mutate(group = factor(cluster,
levels = 1:2,
labels = str_c("n = ", group_size(.)))) %>%
ungroup() %>%
select(participant, cluster, group),
by = "participant")
df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
rating = -10,
group = NA,
participant = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
group = NA,
participant = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = rating,
group = participant,
shape = group)) +
geom_line(mapping = aes(linetype = "individual",
color = group),
size = 1,
alpha = 0.3) +
geom_point(mapping = aes(color = group),
size = 1,
alpha = 0.3) +
stat_summary(mapping = aes(group = group,
color = group,
linetype = "mean"),
fun = "mean",
geom = "line",
size = 1.5) +
stat_summary(data = df.plot %>% filter(cluster == 1),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
stat_summary(data = df.plot %>% filter(cluster == 2),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
annotation_custom(grob = ball_a,
xmin = 0.5,
xmax = 1.5,
ymin = -30,
ymax = -10) +
annotation_custom(grob = ball_b,
xmin = 1.5,
xmax = 2.5,
ymin = -30,
ymax = -10) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_wrap(~ clip,
ncol = 8) +
coord_cartesian(xlim = c(0.9, 2.1),
ylim = c(-5, 105),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25),
labels = seq(0, 100, 25)) +
scale_color_brewer(palette = "Set1",
guide = "none") +
labs(y = "causal responsibility rating",
linetype = "legend",
shape = "") +
theme_bw() +
guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
keywidth = unit(1.2, "cm")),
shape = guide_legend(override.aes = list(color = c("#377eb8",
"#e41a1c"),
shape = c(19, 19),
alpha = c(1, 1),
size = c(5, 5)))) +
theme(panel.grid = element_blank(),
text = element_text(size = 20),
legend.position = c(0.505, 0.23),
axis.title.x = element_blank(),
legend.text = element_text(size = 20),
legend.title = element_text(size = 24),
legend.background = element_rect(fill = "transparent"),
strip.text = element_text(size = 30),
axis.text.x = element_blank(),
legend.key = element_rect(fill = "transparent"),
legend.box = "horizontal",
legend.spacing.x = unit(0.1, "cm"),
panel.spacing.x = unit(1, "cm"),
plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))
print(p)
# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
# plot = p,
# width = 20,
# height = 8.5,
# device = cairo_pdf)
4.5.3 Tables
4.5.3.1 Relationship between predictors
df.exp3.model %>%
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
pivot_longer(cols = A_difference:B_robust,
names_to = c("ball", "aspect"),
names_sep = "_",
values_to = "value") %>%
pivot_wider(names_from = aspect,
values_from = value) %>%
select(difference, whether, how, sufficient, robust) %>%
correlate() %>%
shave() %>%
fashion() %>%
print_table()| rowname | difference | whether | how | sufficient | robust |
|---|---|---|---|---|---|
| difference | |||||
| whether | .50 | ||||
| how | .79 | .27 | |||
| sufficient | .21 | .10 | .36 | ||
| robust | .43 | .93 | .24 | .20 |
4.5.3.2 Population level predictors in the mixed effects models
fit.brm_exp3_causal_w %>%
tidy(effects = "fixed") %>%
mutate(model = "CSM_w") %>%
bind_rows(fit.brm_exp3_causal_wh %>%
tidy(effects = "fixed") %>%
mutate(model = "CSM_wh")) %>%
bind_rows(fit.brm_exp3_causal_whs %>%
tidy(effects = "fixed") %>%
mutate(model = "CSM_whs")) %>%
mutate(term = tolower(term)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
select(model, everything(), -effect, -component) %>%
print_table()| model | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|
| CSM_w | (intercept) | 31.68 | 1.61 | 28.41 | 34.87 |
| CSM_w | whether | 39.30 | 2.16 | 35.11 | 43.54 |
| CSM_wh | (intercept) | 10.29 | 1.56 | 7.30 | 13.41 |
| CSM_wh | whether | 28.40 | 2.68 | 23.18 | 33.56 |
| CSM_wh | how | 36.07 | 2.60 | 31.01 | 41.17 |
| CSM_whs | (intercept) | 10.43 | 1.59 | 7.30 | 13.50 |
| CSM_whs | whether | 27.93 | 2.60 | 22.80 | 33.09 |
| CSM_whs | how | 26.47 | 2.91 | 20.83 | 32.42 |
| CSM_whs | sufficient | 30.90 | 3.06 | 25.03 | 36.93 |
4.5.3.3 Individual participants regression results
df.exp3.causal.individual_fit %>%
select_if(~ (is.numeric(.) | is.factor(.))) %>%
mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>%
select(participant, everything(), best_model) %>%
print_table()| participant | r_w | r_wh | r_whs | rmse_w | rmse_wh | rmse_whs | best_model |
|---|---|---|---|---|---|---|---|
| 1 | 0.29 | 0.38 | 0.48 | 30.00 | 29.03 | 27.74 | whs |
| 2 | 0.20 | 0.77 | 0.77 | 39.32 | 27.94 | 27.55 | whs |
| 3 | 0.37 | 0.78 | 0.79 | 38.90 | 28.00 | 26.95 | whs |
| 4 | 0.41 | 0.55 | 0.63 | 43.62 | 40.60 | 38.31 | whs |
| 5 | 0.45 | 0.53 | 0.54 | 39.23 | 37.32 | 36.81 | whs |
| 6 | 0.18 | 0.54 | 0.54 | 40.84 | 36.32 | 36.07 | whs |
| 7 | 0.49 | 0.61 | 0.64 | 32.96 | 29.93 | 29.15 | whs |
| 8 | 0.39 | 0.63 | 0.68 | 38.37 | 33.26 | 31.64 | whs |
| 9 | 0.40 | 0.72 | 0.76 | 35.86 | 28.34 | 26.39 | whs |
| 10 | 0.28 | 0.52 | 0.53 | 29.13 | 26.28 | 25.85 | whs |
| 11 | 0.51 | 0.56 | 0.60 | 33.66 | 32.19 | 31.17 | whs |
| 12 | 0.56 | 0.57 | 0.70 | 32.78 | 32.08 | 28.84 | whs |
| 13 | 0.56 | 0.71 | 0.74 | 29.83 | 25.40 | 24.21 | whs |
| 14 | 0.43 | 0.47 | 0.50 | 33.72 | 32.82 | 32.23 | whs |
| 15 | 0.19 | 0.35 | 0.37 | 39.77 | 38.29 | 37.94 | whs |
| 16 | 0.62 | 0.67 | 0.76 | 40.49 | 37.98 | 34.39 | whs |
| 17 | 0.42 | 0.66 | 0.71 | 28.57 | 23.82 | 22.42 | whs |
| 18 | 0.31 | 0.35 | 0.39 | 37.78 | 37.09 | 36.57 | whs |
| 19 | 0.37 | 0.60 | 0.62 | 33.41 | 29.23 | 28.62 | whs |
| 20 | 0.57 | 0.61 | 0.68 | 36.47 | 35.00 | 32.71 | whs |
| 21 | 0.43 | 0.50 | 0.57 | 30.95 | 29.61 | 28.32 | whs |
| 22 | 0.51 | 0.67 | 0.71 | 33.72 | 29.36 | 28.00 | whs |
| 23 | 0.33 | 0.45 | 0.52 | 38.76 | 37.00 | 35.57 | whs |
| 24 | 0.46 | 0.63 | 0.69 | 34.01 | 29.96 | 28.22 | whs |
| 25 | 0.61 | 0.66 | 0.70 | 30.60 | 28.91 | 27.55 | whs |
| 26 | 0.57 | 0.71 | 0.76 | 40.43 | 35.03 | 32.80 | whs |
| 27 | 0.46 | 0.52 | 0.66 | 43.12 | 41.49 | 38.12 | whs |
| 28 | 0.30 | 0.45 | 0.46 | 29.33 | 27.62 | 27.38 | whs |
| 29 | 0.48 | 0.60 | 0.65 | 38.93 | 36.03 | 34.48 | whs |
| 30 | 0.18 | 0.70 | 0.70 | 33.06 | 24.95 | 24.71 | whs |
| 31 | 0.48 | 0.59 | 0.69 | 37.54 | 34.72 | 31.76 | whs |
| 32 | 0.37 | 0.61 | 0.70 | 40.62 | 35.94 | 32.92 | whs |
| 33 | 0.50 | 0.56 | 0.64 | 38.64 | 36.79 | 34.85 | whs |
| 34 | 0.24 | 0.51 | 0.52 | 44.50 | 40.86 | 40.19 | whs |
| 35 | 0.15 | 0.32 | 0.30 | 34.56 | 33.36 | 33.37 | wh |
| 36 | 0.19 | 0.84 | 0.82 | 45.38 | 29.58 | 29.70 | wh |
| 37 | 0.55 | 0.63 | 0.67 | 36.04 | 33.36 | 32.11 | whs |
| 38 | 0.52 | 0.61 | 0.67 | 36.03 | 33.39 | 31.53 | whs |
| 39 | 0.47 | 0.77 | 0.77 | 29.61 | 21.75 | 21.58 | whs |
| 40 | 0.46 | 0.78 | 0.79 | 32.41 | 23.60 | 22.69 | whs |
| 41 | 0.32 | 0.50 | 0.58 | 32.90 | 30.38 | 28.69 | whs |
4.5.3.4 CSMwhs predictions for selection of cases
df.exp3.model.causal %>%
left_join(df.exp3.causal.regression,
by = c("clip", "ball")) %>%
filter(clip %in% c(7, 23, 16, 3, 15)) %>%
mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic),
.funs = ~ round(., 2)) %>%
select(clip, ball, difference, whether, how, sufficient, robust) %>%
mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
.funs = ~ ifelse(. < 0.5,
str_c("xmark (", ., ")"),
str_c("cmark (", ., ")"))) %>%
arrange(clip) %>%
print_table()| clip | ball | difference | whether | how | sufficient | robust |
|---|---|---|---|---|---|---|
| 7 | A | cmark (1) | xmark (0.34) | cmark (1) | xmark (0) | xmark (0.25) |
| 7 | B | cmark (1) | cmark (1) | cmark (1) | cmark (0.67) | cmark (0.6) |
| 23 | A | xmark (0.05) | xmark (0) | xmark (0) | xmark (0) | xmark (0) |
| 23 | B | cmark (0.91) | cmark (0.79) | xmark (0) | xmark (0) | cmark (0.72) |
| 16 | A | cmark (1) | xmark (0.23) | cmark (1) | cmark (1) | xmark (0.35) |
| 16 | B | xmark (0) | xmark (0) | xmark (0) | xmark (0) | xmark (0) |
| 3 | A | cmark (1) | cmark (0.88) | cmark (1) | xmark (0.12) | cmark (0.76) |
| 3 | B | cmark (1) | cmark (0.89) | cmark (1) | xmark (0.11) | cmark (0.75) |
| 15 | A | cmark (1) | xmark (0.01) | cmark (1) | cmark (0.99) | xmark (0.1) |
| 15 | B | cmark (1) | xmark (0.01) | cmark (1) | cmark (1) | xmark (0.1) |
4.5.3.5 CSMwhs predictions for all cases
df.exp3.causal.regression %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")) %>%
left_join(df.exp3.clipinfo %>%
select(-clipindex),
by = c("clip")) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
.funs = ~ . * 100) %>%
mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
select(clip, ball, contains("outcome"), difference, whether, how, sufficient,
robust, w, wh, whs, heuristic, rating = rating_mean) %>%
# print_table(format = "latex", digits = 0)
print_table(digits = 0)| clip | ball | outcome_both | outcome_a | outcome_b | outcome_none | difference | whether | how | sufficient | robust | w | wh | whs | heuristic | rating |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | A | 0 | 0 | 0 | 0 | 100 | 40 | 100 | 23 | 36 | 47 | 58 | 55 | 57 | 42 |
| 1 | B | 0 | 0 | 0 | 0 | 100 | 15 | 100 | 16 | 9 | 38 | 51 | 46 | 54 | 37 |
| 2 | A | 0 | 0 | 0 | 0 | 57 | 12 | 0 | 0 | 10 | 36 | 14 | 14 | 25 | 21 |
| 2 | B | 0 | 0 | 0 | 0 | 18 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 24 | 19 |
| 3 | A | 1 | 0 | 0 | 0 | 100 | 88 | 100 | 12 | 76 | 66 | 71 | 65 | 72 | 76 |
| 3 | B | 1 | 0 | 0 | 0 | 100 | 89 | 100 | 11 | 75 | 67 | 72 | 65 | 72 | 75 |
| 4 | A | 1 | 0 | 0 | 0 | 100 | 78 | 100 | 4 | 78 | 62 | 69 | 60 | 58 | 63 |
| 4 | B | 1 | 0 | 0 | 0 | 100 | 95 | 100 | 15 | 57 | 69 | 73 | 68 | 54 | 78 |
| 5 | A | 0 | 0 | 1 | 0 | 100 | 90 | 100 | 0 | 47 | 67 | 72 | 62 | 47 | 71 |
| 5 | B | 0 | 0 | 1 | 0 | 100 | 0 | 100 | 0 | 0 | 32 | 46 | 37 | 68 | 22 |
| 6 | A | 0 | 0 | 1 | 0 | 100 | 59 | 100 | 16 | 35 | 55 | 63 | 59 | 53 | 73 |
| 6 | B | 0 | 0 | 1 | 0 | 100 | 18 | 100 | 6 | 14 | 39 | 52 | 44 | 53 | 22 |
| 7 | A | 1 | 0 | 1 | 0 | 100 | 34 | 100 | 0 | 25 | 45 | 56 | 46 | 70 | 59 |
| 7 | B | 1 | 0 | 1 | 0 | 100 | 100 | 100 | 67 | 60 | 71 | 75 | 86 | 64 | 79 |
| 8 | A | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 25 | 7 |
| 8 | B | 1 | 0 | 1 | 0 | 100 | 100 | 100 | 100 | 100 | 71 | 75 | 96 | 84 | 92 |
| 9 | A | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 14 | 8 |
| 9 | B | 0 | 1 | 0 | 0 | 100 | 100 | 100 | 0 | 100 | 71 | 75 | 65 | 78 | 90 |
| 10 | A | 0 | 1 | 0 | 0 | 77 | 18 | 0 | 0 | 22 | 39 | 15 | 16 | 15 | 23 |
| 10 | B | 0 | 1 | 0 | 0 | 98 | 79 | 0 | 0 | 63 | 63 | 33 | 33 | 21 | 55 |
| 11 | A | 1 | 1 | 0 | 0 | 100 | 70 | 100 | 77 | 68 | 59 | 66 | 80 | 71 | 93 |
| 11 | B | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 16 | 4 |
| 12 | A | 1 | 1 | 0 | 0 | 100 | 82 | 100 | 74 | 83 | 64 | 70 | 83 | 57 | 77 |
| 12 | B | 1 | 1 | 0 | 0 | 100 | 0 | 100 | 12 | 24 | 32 | 46 | 41 | 53 | 37 |
| 13 | A | 0 | 1 | 1 | 0 | 67 | 34 | 0 | 0 | 35 | 45 | 20 | 20 | 14 | 8 |
| 13 | B | 0 | 1 | 1 | 0 | 70 | 35 | 0 | 0 | 35 | 46 | 20 | 20 | 21 | 64 |
| 14 | A | 0 | 1 | 1 | 0 | 97 | 91 | 0 | 0 | 59 | 67 | 36 | 36 | 21 | 22 |
| 14 | B | 0 | 1 | 1 | 0 | 91 | 77 | 0 | 0 | 51 | 62 | 32 | 32 | 20 | 18 |
| 15 | A | 1 | 1 | 1 | 0 | 100 | 1 | 100 | 99 | 10 | 32 | 47 | 68 | 71 | 76 |
| 15 | B | 1 | 1 | 1 | 0 | 100 | 1 | 100 | 100 | 10 | 32 | 47 | 68 | 71 | 76 |
| 16 | A | 1 | 1 | 1 | 0 | 100 | 23 | 100 | 100 | 35 | 41 | 53 | 74 | 80 | 92 |
| 16 | B | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 14 | 4 |
| 17 | A | 0 | 0 | 0 | 1 | 100 | 19 | 100 | 37 | 18 | 39 | 52 | 54 | 66 | 69 |
| 17 | B | 0 | 0 | 0 | 1 | 100 | 0 | 100 | 36 | 17 | 32 | 46 | 48 | 65 | 46 |
| 18 | A | 0 | 0 | 0 | 1 | 100 | 11 | 100 | 40 | 17 | 36 | 49 | 52 | 55 | 63 |
| 18 | B | 0 | 0 | 0 | 1 | 100 | 7 | 100 | 37 | 9 | 34 | 48 | 50 | 56 | 66 |
| 19 | A | 1 | 0 | 0 | 1 | 100 | 74 | 100 | 7 | 65 | 61 | 67 | 60 | 55 | 53 |
| 19 | B | 1 | 0 | 0 | 1 | 100 | 72 | 100 | 7 | 65 | 60 | 67 | 59 | 55 | 49 |
| 20 | A | 1 | 0 | 0 | 1 | 100 | 92 | 100 | 8 | 72 | 68 | 72 | 65 | 57 | 41 |
| 20 | B | 1 | 0 | 0 | 1 | 100 | 88 | 100 | 4 | 53 | 66 | 71 | 63 | 56 | 71 |
| 21 | A | 0 | 0 | 1 | 1 | 100 | 47 | 100 | 40 | 45 | 50 | 60 | 63 | 58 | 80 |
| 21 | B | 0 | 0 | 1 | 1 | 100 | 9 | 100 | 21 | 10 | 35 | 49 | 46 | 59 | 18 |
| 22 | A | 0 | 0 | 1 | 1 | 100 | 100 | 100 | 89 | 83 | 71 | 75 | 92 | 48 | 60 |
| 22 | B | 0 | 0 | 1 | 1 | 100 | 8 | 100 | 0 | 15 | 35 | 49 | 39 | 53 | 42 |
| 23 | A | 1 | 0 | 1 | 1 | 5 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 15 | 3 |
| 23 | B | 1 | 0 | 1 | 1 | 91 | 79 | 0 | 0 | 72 | 63 | 33 | 33 | 22 | 39 |
| 24 | A | 1 | 0 | 1 | 1 | 100 | 66 | 100 | 4 | 63 | 58 | 65 | 57 | 57 | 44 |
| 24 | B | 1 | 0 | 1 | 1 | 100 | 94 | 100 | 22 | 79 | 69 | 73 | 70 | 54 | 73 |
| 25 | A | 0 | 1 | 0 | 1 | 100 | 25 | 100 | 21 | 26 | 42 | 53 | 50 | 69 | 43 |
| 25 | B | 0 | 1 | 0 | 1 | 100 | 74 | 100 | 54 | 65 | 61 | 67 | 74 | 56 | 73 |
| 26 | A | 0 | 1 | 0 | 1 | 100 | 6 | 100 | 3 | 9 | 34 | 48 | 40 | 60 | 39 |
| 26 | B | 0 | 1 | 0 | 1 | 100 | 87 | 100 | 35 | 54 | 66 | 71 | 72 | 46 | 69 |
| 27 | A | 1 | 1 | 0 | 1 | 100 | 97 | 100 | 52 | 97 | 70 | 74 | 80 | 67 | 80 |
| 27 | B | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 17 | 6 |
| 28 | A | 1 | 1 | 0 | 1 | 100 | 90 | 100 | 22 | 80 | 67 | 72 | 69 | 79 | 89 |
| 28 | B | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 12 | 5 |
| 29 | A | 0 | 1 | 1 | 1 | 100 | 58 | 100 | 24 | 44 | 55 | 63 | 61 | 66 | 47 |
| 29 | B | 0 | 1 | 1 | 1 | 100 | 63 | 100 | 24 | 38 | 57 | 64 | 62 | 54 | 67 |
| 30 | A | 0 | 1 | 1 | 1 | 100 | 57 | 100 | 29 | 49 | 54 | 63 | 62 | 63 | 58 |
| 30 | B | 0 | 1 | 1 | 1 | 100 | 46 | 100 | 24 | 39 | 50 | 60 | 57 | 63 | 56 |
| 31 | A | 1 | 1 | 1 | 1 | 100 | 2 | 100 | 4 | 3 | 32 | 47 | 39 | 63 | 44 |
| 31 | B | 1 | 1 | 1 | 1 | 100 | 4 | 100 | 4 | 4 | 33 | 47 | 39 | 53 | 46 |
| 32 | A | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 16 | 5 |
| 32 | B | 1 | 1 | 1 | 1 | 100 | 75 | 100 | 66 | 73 | 61 | 68 | 78 | 65 | 71 |
4.5.3.6 Heuristic model
fit.brm_exp3_causal_heuristic %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | 49.79 | 1.46 | 46.91 | 52.74 |
| fixed | cond | moving | 0.22 | 0.22 | 0.01 | 0.81 |
| fixed | cond | speed | 2.06 | 0.85 | 0.44 | 3.67 |
| fixed | cond | contact_e | 0.38 | 0.35 | 0.01 | 1.29 |
| fixed | cond | e_speed_diff | 0.12 | 0.12 | 0.00 | 0.46 |
| fixed | cond | e_direction_diff | 1.04 | 0.74 | 0.04 | 2.71 |
| fixed | cond | total_speed_diff | 2.21 | 0.95 | 0.41 | 4.17 |
| fixed | cond | total_direction_diff | 3.99 | 0.91 | 2.14 | 5.68 |
| fixed | cond | transfer | 15.59 | 0.79 | 14.06 | 17.14 |
| fixed | cond | e_moving | 0.18 | 0.17 | 0.01 | 0.62 |
| fixed | cond | exclusive | 4.39 | 0.73 | 2.97 | 5.81 |
5 Appendix
5.1 Experiment 3
5.1.1 Does the question framing matter? (Responsibility vs. causation)
We have run a pilot eye-tracking experiment with 15 participants who saw the same clips as the ones in Experiment 3. Instead of asking participants to judge “To what extent were A and B responsible for E (not) going through the gate?”, we asked them to judge to what extent they agreed with the following statement: “Ball A/B caused ball E to go through the gate.” or “Ball A/B prevented ball E from going through the gate.” depending on the outcome.
Participants’ agreement judgments with the causal statements in this pilot eye-tracking experiment were strikingly similar to participants’ responsibility judgments in the online experiment.
df.exp3.eye_tracking_pilot =
read_csv("../../data/experiment3_eye_tracking_pilot.csv") %>%
rename(clip = trial) %>%
mutate(ball = str_to_upper(ball)) %>%
group_by(clip, ball, outcome) %>%
summarize(causality_mean = mean(causation),
causality_low = smean.cl.boot(causation)[2],
causality_high = smean.cl.boot(causation)[3]) %>%
ungroup()
df.plot = df.exp3.eye_tracking_pilot %>%
# select(clip, ball, outcome,) %>%
left_join(df.exp3.causal.means %>%
select(clip,
ball,
responsibility_mean = rating_mean,
responsibility_low = rating_low,
responsibility_high = rating_high),
by = c("clip", "ball"))
# highlight the double prevention clip (#23)
df.plot = df.plot %>%
mutate(color = ifelse(clip == 23, "1", "0"))
ggplot(data = df.plot,
mapping = aes(x = causality_mean,
y = responsibility_mean)) +
geom_abline(intercept = 0,
slope = 1,
linetype = 2) +
geom_linerange(mapping = aes(ymin = responsibility_low,
ymax = responsibility_high),
alpha = 0.1) +
geom_linerange(mapping = aes(xmin = causality_low,
xmax = causality_high),
alpha = 0.1) +
geom_point(size = 2,
mapping = aes(color = color),
show.legend = F) +
geom_smooth(method = "lm",
color = "black") +
annotate(geom = "text",
x = c(0, 0),
y = c(100, 90),
label = c(str_c("r = ",
cor(df.plot$causality_mean,
df.plot$responsibility_mean) %>%
round(2)),
str_c("RMSE = ", rmse(df.plot$causality_mean,
df.plot$responsibility_mean) %>%
round(2))),
hjust = 0,
size = 8) +
scale_x_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_y_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_color_manual(values = c("black", "red")) +
labs(x = "causal judgment",
y = "responsibility judgment") +
theme(text = element_text(size = 24))
# ggsave("../../figures/plots/aux_question_framing_comparison.pdf",
# width = 6,
# height = 6)
5.1.2 CSM_WHS fits for each individual participant
df.plot = fit.brm_exp3_causal_whs %>%
fitted() %>%
as_tibble() %>%
select(prediction = Estimate) %>%
bind_cols(df.exp3.causal.long) %>%
relocate(prediction, .after = last_col())
df.text = df.plot %>%
group_by(participant) %>%
summarize(r = round(cor(prediction, rating), 2)) %>%
ungroup() %>%
mutate(r = str_c("r = ", r),
prediction = 1,
rating = 110)
ggplot(data = df.plot,
mapping = aes(x = prediction,
y = rating)) +
geom_abline(intercept = 0,
slope = 1,
color = "blue",
alpha = 0.5) +
geom_point(alpha = 0.3,
size = 1) +
geom_text(data = df.text,
mapping = aes(label = r),
hjust = 0,
color = "red",
size = 4) +
facet_wrap(facets = vars(participant),
ncol = 7) +
labs(x = "model prediction",
y = "causal responsibility rating") +
coord_cartesian(xlim = c(0, 100),
ylim = c(0, 100),
expand = F,
clip = "off") +
scale_x_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_y_continuous(breaks = seq(0, 100, 25),
expand = expansion(mult = c(0, 0))) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_blank(),
text = element_text(size = 16),
panel.spacing.x = unit(0.75, "cm"),
panel.spacing.y = unit(1, "cm"),
plot.margin = margin(t = 0.7,
r = 0.8,
b = 0.2,
l = 0.2,
unit = "cm"))
# ggsave("../../figures/plots/exp3_individual_scatter_plots.pdf",
# width = 12,
# height = 8)
5.1.3 Individual participant fit model comparison
df.plot = df.exp3.causal.individual_fit %>%
select(participant, contains("r_"), contains("rmse_")) %>%
pivot_longer(cols = -participant,
names_to = c("measure", "model"),
names_sep = "_",
values_to = "fit")
ggplot(data = df.plot %>%
filter(measure == "r"),
mapping = aes(x = fit,
fill = model)) +
geom_density(alpha = 0.5) +
labs(x = "correlation value") +
scale_x_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
scale_fill_brewer(palette = "Set1") +
# ggplot2::theme_classic() +
theme(legend.position = c(0.3, 0.95),
legend.direction = "horizontal",
text = element_text(size = 20))
# ggsave("../../figures/plots/exp3_individual_densities.pdf",
# width = 8,
# height = 6)
### 3D scatter plot of participant clusters
plot_ly(x = df.cluster_whs$whether,
y = df.cluster_whs$how,
z = df.cluster_whs$sufficient,
type = "scatter3d",
mode = "markers",
color = as.factor(df.cluster_whs$cluster),
colors = c("#e41a1c", "#377eb8")) %>%
layout(showlegend = FALSE,
title = "Participant clusters",
scene = list(
xaxis = list(title = "whether"),
yaxis = list(title = "how"),
zaxis = list(title = "sufficient")))5.1.4 Clusters of participants
set.seed(1)
df.plot = df.exp3.causal.long %>%
left_join(df.cluster_whs %>%
select(participant, cluster) %>%
mutate(participant = as.numeric(participant),
cluster = as.factor(cluster)) %>%
group_by(cluster) %>%
mutate(label = factor(cluster,
levels = 1:2,
labels = str_c("n = ", group_size(.)))) %>%
ungroup(),
by = "participant")
ggplot(data = df.plot,
mapping = aes(x = ball,
y = rating,
group = label,
fill = label)) +
geom_point(shape = 21,
alpha = 0.2,
position = position_jitterdodge(jitter.width = 0.1,
jitter.height = 0,
dodge.width = 0.75)) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
mapping = aes(color = label),
size = 1,
position = position_dodge(width = 0.75)) +
stat_summary(fun = "mean",
geom = "point",
size = 4,
shape = 21,
position = position_dodge(width = 0.75)) +
facet_wrap(facets = vars(clip), ncol = 8) +
labs(y = "causal responsibility rating",
fill = "cluster",
color = "cluster") +
scale_fill_brewer(palette = "Set1") +
scale_color_brewer(palette = "Set1") +
theme(text = element_text(size = 24),
legend.position = "bottom")
# ggsave("../../figures/plots/exp3_clusters_points.pdf",
# width = 12,
# height = 8)
5.1.5 Ternary plots for individual regression results
library("ggtern")
df.plot = df.exp3.causal.individual_fit %>%
select(participant, fit_whs) %>%
mutate(estimates = map(fit_whs, ~ fixef(.) %>%
as_tibble(rownames = "term") %>%
clean_names())) %>%
select(participant, estimates) %>%
unnest(estimates) %>%
filter(!str_detect(term, "Intercept")) %>%
select(participant, term, estimate) %>%
pivot_wider(names_from = term,
values_from = estimate) %>%
mutate_at(vars(whether, how, sufficient),
.funs = list(norm = ~ . / (how + whether + sufficient))) %>%
mutate(color = 0,
color = ifelse(how_norm == max(how_norm), 1, color),
color = ifelse(whether_norm == max(whether_norm), 2, color),
color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
color = factor(color))
ggplot(data = df.plot,
mapping = aes(x = how,
y = sufficient,
z = whether,
color = color)) +
geom_point(alpha = 0.7,
size = 2,
show.legend = F) +
scale_color_manual(values = c("black", "red", "green", "blue")) +
coord_tern() +
theme_showarrows() +
theme(text = element_text(size = 20),
tern.axis.title.T = element_blank(),
tern.axis.title.L = element_blank(),
tern.axis.title.R = element_blank())
# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
# width = 5,
# height = 5)6 Session info
R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.0 purrr_0.3.4
[5] readr_1.3.1 tidyr_1.1.0 tibble_3.0.1 tidyverse_1.3.0
[9] patchwork_1.0.0 clValid_0.6-7 cluster_2.1.0 egg_0.4.5
[13] gridExtra_2.3 plotly_4.9.2.1 png_0.1-7 tidybayes_2.1.1
[17] Hmisc_4.4-0 ggplot2_3.3.1 Formula_1.2-3 survival_3.1-8
[21] lattice_0.20-38 janitor_2.0.1 jsonlite_1.6.1 xtable_1.8-4
[25] corrr_0.4.2 brms_2.13.0 Rcpp_1.0.4.6 broom.mixed_0.2.6
[29] lme4_1.1-23 Matrix_1.2-18 kableExtra_1.1.0 knitr_1.28
loaded via a namespace (and not attached):
[1] readxl_1.3.1 backports_1.1.7 plyr_1.8.6
[4] igraph_1.2.5 lazyeval_0.2.2 svUnit_1.0.3
[7] TMB_1.7.16 splines_3.6.3 crosstalk_1.1.0.1
[10] rstantools_2.1.0 inline_0.3.15 digest_0.6.25
[13] htmltools_0.4.0 rsconnect_0.8.16 fansi_0.4.1
[16] magrittr_1.5 checkmate_2.0.0 modelr_0.1.8
[19] RcppParallel_5.0.1 matrixStats_0.56.0 xts_0.12-0
[22] prettyunits_1.1.1 jpeg_0.1-8.1 colorspace_1.4-1
[25] blob_1.2.1 rvest_0.3.5 ggdist_2.1.1
[28] haven_2.3.1 xfun_0.14 callr_3.4.3
[31] crayon_1.3.4 zoo_1.8-8 glue_1.4.1
[34] gtable_0.3.0 emmeans_1.4.7 webshot_0.5.2
[37] pkgbuild_1.0.8 rstan_2.19.3 abind_1.4-5
[40] scales_1.1.1 mvtnorm_1.1-1 DBI_1.1.0
[43] miniUI_0.1.1.1 viridisLite_0.3.0 htmlTable_1.13.3
[46] HDInterval_0.2.2 foreign_0.8-75 stats4_3.6.3
[49] StanHeaders_2.21.0-5 DT_0.13 htmlwidgets_1.5.1
[52] httr_1.4.1 threejs_0.3.3 arrayhelpers_1.1-0
[55] RColorBrewer_1.1-2 acepack_1.4.1 ellipsis_0.3.1
[58] farver_2.0.3 pkgconfig_2.0.3 loo_2.2.0
[61] dbplyr_1.4.4 nnet_7.3-12 utf8_1.1.4
[64] labeling_0.3 tidyselect_1.1.0 rlang_0.4.6
[67] reshape2_1.4.4 later_1.1.0.1 cellranger_1.1.0
[70] munsell_0.5.0 tools_3.6.3 cli_2.0.2
[73] generics_0.0.2 broom_0.7.0.9000 ggridges_0.5.2
[76] evaluate_0.14 fastmap_1.0.1 yaml_2.2.1
[79] fs_1.4.1 processx_3.4.2 nlme_3.1-144
[82] mime_0.9 xml2_1.3.2 compiler_3.6.3
[85] bayesplot_1.7.2 shinythemes_1.1.2 rstudioapi_0.11
[88] reprex_0.3.0 statmod_1.4.34 stringi_1.4.6
[91] highr_0.8 ps_1.3.3 Brobdingnag_1.2-6
[94] nloptr_1.2.2.1 markdown_1.1 shinyjs_1.1
[97] vctrs_0.3.1 pillar_1.4.4 lifecycle_0.2.0
[100] bridgesampling_1.0-0 estimability_1.3 data.table_1.12.8
[103] httpuv_1.5.4 R6_2.4.1 latticeExtra_0.6-29
[106] bookdown_0.19 promises_1.1.1 boot_1.3-24
[109] colourpicker_1.0 MASS_7.3-51.5 gtools_3.8.2
[112] assertthat_0.2.1 withr_2.2.0 shinystan_2.5.0
[115] mgcv_1.8-31 parallel_3.6.3 hms_0.5.3
[118] rpart_4.1-15 class_7.3-15 coda_0.19-3
[121] minqa_1.2.4 rmarkdown_2.2 snakecase_0.11.0
[124] shiny_1.4.0.2 lubridate_1.7.9 base64enc_0.1-3
[127] dygraphs_1.1.1.6